function fig1_10

%  solve thermokinetic IVP using MATLAB routines

clf

ep=0.1;
tmax=5;

[t,y] = ode45(@rhs,[0 tmax],[0 0]);

% % get(gcf)
% set(gcf,'Position', [7 1140 577 197])

tt=linspace(0,tmax,800);
for it=1:800
    u0(it)=1-exp(-tt(it));
    u1(it)=u0(it)+ep*( 0.5*(tt(it)^2+2*tt(it)-4)*exp(-tt(it)) + (2+tt(it))*exp(-2*tt(it))) ;
    q0(it)=1-(1+tt(it))*exp(-tt(it));
    q1(it)=q0(it)+ep*( (tt(it)^3-18*tt(it)+30)*exp(-tt(it))/6 - (2*tt(it)+5)*exp(-2*tt(it))) ;
end

% figure(1)
% hold on
% box on
% yy=y(:,1);
% plot(t,yy)
% plot(tt,u0,'--r','LineWidth',1)
% plot(tt,u1,'-.g','LineWidth',1)
% legend(' Exact',' 1 term',' 2 terms','Location','SouthEast','FontSize',16)
% xlabel('t-axis','FontWeight','bold')
% ylabel('Solution','FontWeight','bold')
% set(gca,'FontSize',16)
% hold off

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.2,'MR',-0.01,'ML',0.05,'P',0.04)

box on
grid on
yy=y(:,2);
plot(tt,q0,'-.k','LineWidth',1.2)
hold on
plot(tt,q1,'r','LineWidth',1.2)
plot(t,yy,'--b','LineWidth',1.4)
axis([0 tmax 0 1])
ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$t$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('$q(t)$','Interpreter','latex','FontSize',14)

legend(' 1 term',' 2 terms',' Numerical','Location','SouthEast','FontSize',13,'Interpreter','latex')
hold off

%exportgraphics(gcf,'/Users/mark/Desktop/q.pdf')

function dy=rhs(t,y)
dy=zeros(2,1);
ep=0.1;
dy(1)=1-y(1)*exp(ep*(y(2)-1));
dy(2)=y(1)*exp(ep*(y(2)-1)) - y(2);











